UNCLASSIFIED 


Defense  Technical  Information  Center 
Compilation  Part  Notice 

ADPO 13743 

TITLE:  Knot  Removal  for  Tensor  Product  Splines 
DISTRIBUTION:  Approved  for  public  release,  distribution  unlimited 


This  paper  is  part  of  the  following  report: 

TITLE:  Algorithms  For  Approximation  IV.  Proceedings  of  the  2001 
International  Symposium 

To  order  the  complete  compilation  report,  use:  ADA412833 

The  component  part  is  provided  here  to  allow  users  access  to  individually  authored  sections 
of  proceedings,  annals,  symposia,  etc.  However,  the  component  should  be  considered  within 
the  context  of  the  overall  compilation  report  and  not  as  a stand-alone  technical  report. 

The  following  component  part  numbers  comprise  the  compilation  report: 

ADPO  13708  thru  ADPO  13761 


UNCLASSIFIED 


Knot  removal  for  tensor  product  splines 

T.  Brenna 

Dept,  of  Informatics,  Univ.  of  Oslo,  Oslo. 
trondbreSifi.uio.no 


Abstract 

Given  a spline  function  as  a B-spline  expansion  the  object  of  knot  removal  is  to  remove  as 
many  knots  as  possible  without  perturbing  the  spline  by  more  than  a specified  tolerance. 
In  1987  Lyche  and  Mprken  proposed  an  efficient  knot  removal  algorithm  which  determines 
both  the  number  of  remaining  knots  and  their  position  automatically.  In  this  paper  we 
show  how  their  method  can  be  extended  to  knot  removal  techniques  for  multivariate 
tensor  product  splines.  We  propose  a number  of  new  strategies  for  removing  as  many 
knots  as  possible,  and  discuss  some  of  the  advantages  and  challenges  posed  by  the  special 
structure  of  tensor  product  splines. 


1 Introduction 

Given  a spline  function  we  are  often  interested  in  an  approximate  representation  re- 
quiring less  data.  The  object  of  knot  removal  is  to  remove  as  many  knots  as  possible 
from  a given  spline  without  perturbing  the  spline  by  more  than  a given  tolerance.  An 
efficient  knot  removal  strategy  presented  in  [6]  determines  both  the  number  of  remaining 
knots  and  their  location  automatically.  This  strategy  was  later  extended  to  parametric 
curves  and  surfaces  in  [5],  and  incorporated  with  various  constraints  such  as  monoton- 
icity and  convexity  in  [1].  An  efficient  implementation  of  knot  removal  for  the  special 
case  of  trilinear  splines  is  given  in  [3].  In  this  paper  we  address  some  of  the  questions 
and  problems  arising  when  extending  the  knot  removal  technique  to  multivariate  tensor 
product  splines. 

The  outline  of  this  paper  is  as  follows.  We  start  by  fixing  notation  and  presenting 
techniques  for  representing  tensor  product  splines.  We  then  proceed  with  generalizations 
of  coefficient  norms,  approximation  methods,  methods  for  ranking  the  knots  etc.,  as  we 
review  the  central  parts  of  the  knot  removal  strategy.  Two  different  ways  of  performing 
knot  removal  are  given  together  with  accompanying  strategies  for  finding  the  desired 
approximations.  We  end  the  paper  with  two  examples  demonstrating  various  aspects  of 
the  knot  removal  techniques  presented. 

2 Notation 

Let  d = (dk),  m = (mk)  € Zs  with  0 < d < m (component-wise)  for  some  positive 
integer  s.  Also  let  tk  = {t,-}”Dl+dk+}  be  a knot  vector  with  dk  + 1 equal  knots  at  both 
ends  and  with  no  knot  value  occurring  more  than  dk  + 1 times,  for  k = 1, . . .,  s.  In  this 
paper  we  will  treat  the  collection  t = {tfc}£=1  as  a “single”  knot  vector  with  “length” 
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m + d + 1 defined  to  be  the  sum  of  the  length  of  the  knot  vectors  tfc.  k = 1, . . s.  Given 
such  a knot  vector  we  may  form  products  of  the  basis  functions  associated  with  each 
individual  knot  vector  tk.  By  letting 

S 

A(x)  = Bi]d)t(x)  = Bik4k^(xk)  for  1 < i < m, 
fe= i 

where  i = (ik)  € Zs  and  x = (xk)  £ Ks  ,we  get  a total  of  Ofe=i  TOfc  new  basis  functions 

5 

for  the  tensor  product  space  §d,t  = ® &dk,tk ■ In  this  paper  we  let  B,.  dk  tk  be  the  zfcth 

B-spline  of  degree  dk  associated  with  tfc,  for  k = 1, . . s. 

To  represent  an  element  of  §d,t  we  use  a variant  of  the  classical  Kronecker  product  of 
matrices.  Recall  that  if  A = (ay)™1^^  6 B = (by)”21’"i1  e then  this 

product  is  given  by  A ® B = (aijBjj^’jii*  In  this  paper  we  will  use  the  “equivalent” 
product  defined  by  A ® B = (Ab|j)“^’jf.1,  which  gives  a more  convenient  ordering  of 
the  matrix  elements  for  our  use.  Also  recall  that  for  real  matrices  A,B,C,D  we  have  the 
following  useful  relations  (assuming  that  the  matrix  products  and  inverses  are  defined) 
(A®B)(C®D)  = (AC)®  (BD),  (A®B)“1  = A^B"1  and  A®B  = PX(B®  A)P2, 
for  some  permutation  matrices  Pi  and  P2.  In  addition  we  have  that  the  product  A®B 
will  have  linearly  independent  columns,  provided  the  same  holds  for  A and  B.  For  further 
properties  of  the  Kronecker  product  we  refer  to  [4], 

An  element 

XtX  l TTL  g S 

/(*)  = E • • • E n = E e §d,t 

i\—l  is= 1 k= 1 i<m 

can  now  be  written 

/(x)  = B?f, 

where  Bt  = ® Btk  with  Btk  = (Bl  dkitk, . . .,Bmkidk,tk)T  for  k = 1 Here  f is 

k— 1 

a vector  containing  the  B-spline  coefficients  F = (fi1(..,is)  of  / given  by  f = vec(F)  := 

S 

2^i<m  bei,  where  ei  = ® e,k  with  eik  e Wnk . Finally  we  state  that  for  a tensor  of  real 
~ k= 1 

coefficients  F = (fi)i<i<m  € Em  we  let  denote  the  tensor  F with  its  elements 

rearranged  according  to  the  cyclic  permutation  of  the  s-tuple  {1, 2, . . .,  s}  given  by  ak  = 
{fc,  k + 1, . . .,s,  1, . . .,k  — 1},  for  k = 1, . . .,  s. 

Finally,  for  a spline  / = Xu-cm  /i-®i,d,t(x)  we  define  a class  of  weighted  Zp-norms  of 
its  B-spline  coefficients,  given  by 


_/(Ei<mu;il/ilP)1/p- 
l|/lkt-l  mas  I/ll, 

v l<i<m 


for 

for 


1 < p < oo, 
P = oo,' 


t * —t* 

where  the  weights  are  given  by  w\  = JIfc=i  'fc+d*+i  ~ k > for  1 < i < m.  Using  the 
notation  introduced  above  we  have  that  ||  / |!(p,t=||  Wypf  ||iP,  (p  > 1)  where  Wt  is  a 
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diagonal  scaling  matrix  given  by 

W*  = J1W1‘.  with  W„  - diag  ( (%^)  , ■ f ^ ■))  • 

These  coefficient  norms  are  easy  to  compute  and  are  known  to  approximate  the  ordinary 
Lp -norms  well  for  splines  of  moderate  degree  [2,6].  In  the  algorithms  we  use  p = 2 when 
computing  approximations  and  p = oo  to  measure  the  error. 

3 The  knot  removal  algorithm 

Given  an  element  / G §d,t,  a tolerance  e > 0 and  some  norm  ||  • ||  the  goal  of  the  knot 
removal  algorithm  presented  in  [6]  is  to  find  a subspace  §d.T  of  §d,t  (t  C t)  and  an 
element  g G §d,r  with  ||  / — g ||  < e,  and  where  we  want,  r to  be  of  minimal  length.  In 
this  section  we  review  the  basic  parts  of  this  algorithm  as  we  extend  the  theory  to  tensor 
product  splines.  Further  details  of  the  material  in  this  section  can  be  found  in  [2]. 

3.1  Finding  approximations 

To  approximate  / G §d,t  in  a subspace  §d,T>  where  r is  of  “length”  n + d + 1 with 
n < m,  we  use  the  spline  g which  is  the  best  approximation  to  / in  the  /2,  t-norm. 
In  other  words,  the  spline  we  seek  will  be  the  solution  to  the  minimization  problem 
^min  ||  / - h ||22  t.  Solving  this  problem  is  equivalent  to  solving  the  linear  least  squares 

problem  given  by 

min  ||wt1/2(Ac-f)||?2,  (3.1) 

S 

where  A = ® Ak  is  the  knot  insertion  matrix  from  r to  t (i.e.  Ak  is  the  knot  insertion 
k=l 

matrix  from  Tk  to  tfc,  for  k = 1, . . .,s),  f = vec(F)  are  the  given  B-spline  coefficients 
of  / in  §d,t  and  c = vec(C)  are  the  unknown  B-spline  coefficients  of  g in  §d.T-  Since 
the  knot  insertion  matrix  A has  full  rank  and  Wt  is  non-singular,  the  normal  equations 
ATWtAc  = ATWtf  associated  with  the  system  (3.1)  will  have  a unique  solution  which 
can  be  found  ([2,3])  by  solving  a series  of  s tensor  equation  systems  given  by 

(AjWtkAk)D^5  = (A^  Wtk)Djc'^kj , (3.2) 

for  k = 1, . . s.  Here  Dk  G M”k  with  n*.  = (ni, . . 1, . . ms),  and  we  let  Do  = 

F,  and  set  the  coefficients  of  the  approximation  g equal  to  the  solution  of  the  last 
tensor  equation  system,  C = Ds.  The  tensor  equations  (3.2)  can  be  efficiently  solved 
by  calculating  the  Cholesky  factorization  of  the  banded  coefficient  matrix  (A£WtkAk) 
and  solving  for  each  right  hand  side  in  the  tensor  (A^W^jD^j. 

3.2  Ranking  the  knots 

The  final  approximation  to  the  initial  spline  is  found  by  searching  through  a sequence  of 
approximations,  constructed  by  using  the  approximation  method  of  the  previous  section, 
on  subsets  of  the  knots  of  the  initial  spline.  These  subsets  are  calculated  by  associating  a 
weight  with  each  interior  knot,  representing  a rough  measure  of  its  importance.  See  [6]  for 
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the  details.  For  higher  dimensional  tensor  product  splines  we  set  the  weight  for  a given 
knot  to  the  maximum  of  the  weights  corresponding  to  this  knot  when  the  calculation  is 
iterated  over  the  “remaining”  parameter  directions.  We  refer  to  [2]  for  further  details. 

4 Knot  removal  methods 


When  removing  knots  from  a tensor  product  spline  we  are  faced  with  more  options  than 
in  the  case  of  a spline  curve.  In  this  section  we  present  two  different  ways  of  performing 
knot  removal.  The  first  one  studied  in  [2]  based  on  a symmetric  approach,  treats  all  the 
parameter  directions  of  a tensor  product  spline  simultaneously,  while  the  second  one  will 
treat  one  parameter  direction  at  a time. 

4.1  Knot  removal  based  on  a symmetric  approach 

If  we  let  Gf(r)  denote  the  approximation  to  / G §d,t  defined  on  the  knot  vector 
r we  see  that  the  approximations  in  the  sequence  mentioned  above  can  be  written 
{G fir where  Tj  is  constructed  from  t by  removing  j of  its  interior  knots,  and 
N = iimk  ~ (dk  + 1)]  is  the  total  number  of  interior  knots  of  t.  Given  such  a 
sequence  of  approximations  we  can  perform  a search  on  the  index  j to  determine  an 
approximation  g*  = Gf(r*)  to  the  initial  spline  / with  a preferably  short  knot  vector 
t*,  and  with  the  property  that  ||  f — g*  ||z°°,t<  e~,  where  e is  the  specified  tolerance.  If  the 
knot  vector  t*  is  not  equal  to  any  of  the  two  knot  vectors  r0  or  t/v  we  may  repeat  the 
process  to  find  a new  approximation  based  on  g*  as  proposed  in  [6].  Taking  into  account 
how  the  sequence  {Gf(rj)}jL0  was  constructed  we  expect  the  error  ||  / — Gf(rj)  ||;°o)t 
to  decrease,  but  not  necessarily  strictly,  for  decreasing  values  of  the  search  parameter 
j.  How  the  search  among  the  possible  approximations  is  done  will  generally  depend  on 
a number  of  factors,  including  some  which  will  be  discussed  later  through  examples. 
Also  note  that  we  only  have  to  compute  approximations  for  indexes  actually  used  in 
the  search.  By  treating  all  the  directions  simultaneously  we  take  into  consideration  the 
inherent  symmetry  of  the  problem.  As  we  will  see  later  this  will  in  some  cases  enable 
us  to  remove  more  knots  than  by  treating  one  parameter  direction  at  a time,  but  it  will 
also  lead  to  more  complicated  and  slower  code  in  an  implementation. 


4.2  Knot  removal  for  one  parameter  direction  at  a time 


In  the  second  knot  removal  method  we  start  by  thinking  of  a spline  / € §d,t  as  a series 
of  parametric  curves  in  corresponding  high  dimensional  spaces.  We  can  then  perform  a 
parametric  knot  removal  for  each  parameter  direction.  The  advantage  of  this  approach 
is  that  it  is  easy  to  implement  since  we  may  use  existing  knot  removal  routines  for  spline 
curves  with  only  minor  modifications. 

In  the  following  discussion  we  let  e — X^=i  G,  with  e*  > 0 for  all  i,  be  a given 


tolerance.  Also  let  /(x)  = ^Ti<m  /iBi,d,t(x)  = B^f  be  a spline  in  §d,t 


® § 
fc=i 


dk, tk  j 


with 


B?  = ,®BJ  and  f = vec(F).  We  start  by  identifying  a series  of  parametric  curves 
k=l 

which  may  be  naturally  associated  with  this  tensor  product  spline.  We  say  that  the 
spline  / consists  of  the  curves  fk{xk),  for  k = l,...,s,  where  fk(xk)  is  the  parametric 
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curve  in  RMk,  for  Mk  = (n^=i  ^uOdl^r+i  m*)>  Siven  b>r 
fk{xk)  = ® Bjl  ® (j  ® 


We  now  return  to  the  problem  of  finding  a preferably  short  knot  vector  rCt  and  a spline 
*(*)  = Sj<n  cjf?j,d,r(x)  e SdiT  = ®Sdk  Tk  with  the  property  that  ||  / - g ||j°°,t<  e.  To 
apply  knot  removal  to  / G §djt  we  can  now  go  through  the  following  steps  for  k = 1, . . s. 


1.  Apply  parametric  knot  removal  with  the  tolerance  ek  to  the  parametric  curve 

fk(xk)  = ( ® Ini')  ® ® ( ® Imi^l  fk— 1> 

L\/=l  / \l=k+l  /J 

defined  on  tk , starting  with  f0  = f . 

2.  This  will  produce  a new  parametric  curve  defined  on  the  knot  vector  rk  C tk 

fk(xk)  — [((®  Ini)  ®B^k  ® 

where  fk  — vec(Fk)  for  Fk  € K*n>->“k,n»k+ii-in»»# 

3.  We  also  have  that 


fk, 


= 

r/k-i  \ 

( ® Ini  ] 

I ® B Jk  ® ( 

® Imi^jfk 

L M=1  ) 

U=k+1  /J 

' fk- 1 \ 

| ® b£  ® ( 

.l=l+lImi)]  [(,§/"■)  ( 

® Im,)]fk 
U=k+1  / 

where  Ak  is  the  knot  insertion  matrix  from  rfc  to  tk . 
4.  And  consequently 

/fc-i 


II  fk  ~ f'k  lli«\t't  — 


ffc-1  - [(  ®jln,)  ® Ak  ® ^_® 


fk 


< ek. 


Finally  we  let  the  coefficients  of  the  function  g(x)  = B^c  e §d,T  be  c = vec(Fs),  and 
we  have  the  following  result. 


Theorem  4.1  If  we  let  /(x)  = BJFT  e §d;t  and  g(x)  = B^c  € §d,T  be  the  tensor 
product  splines  from  the  discussion  above,  then  we  have  ||  / — g ||/“,t<  £• 


Proof:  Let  A = ® Ak  be  the  knot  insertion  matrix  from  r to  t,  and  let  /o(x)  = B^f0 
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be  equal  to  / and  /s(x)  = B^fs  be  equal  to  g,  i.e.  fo  = f and  fs  = c.  Then 


\f~9\ 


\t=  ||  IQ 


fn  - Afs 


\l°° 


fo  + XX[(  ;?1Al)  0 - [(«>»  ® (jfjm^fk-l)  - ® Akfs 

k=2  ~ ~ 

s Ell  [(*>.)«(,!  Imi^j  [fk-1  0 -^k  0 (i-0-fi^mi)] 


k= 1 


fk 


ffc-l  _ (j0lni)  0 Ak®  (j_®  ^m^jfk  =Xl||  fk~fi 

k=  1 ~ 1 — +1  ;oo  fe=1 


<£.  □ 


5 Examples 

The  knot  removal  methods  presented  above  have  been  implemented  and  tested  on  a 
computer.  In  this  section  we  present  trivariate  examples  from  this  implementation  and 
propose  different  knot  removal  strategies  depending  on  the  problem  at  hand.  See  [3]  for 
a detailed  description  of  this  implementation. 

Example  5.1  In  this  first  example  we  will  compare  two  different  strategies  for  searching 
through  a list  of  approximations  {G/(rj)}^L 0 introduced  above.  We  will  consider  the 
knot  removal  method  treating  one  parameter  direction  at  a time,  which  means  that  we 
end  up  solving  a parametric  knot  removal  problem  with  tolerance  gj  = e/3,  * = 1,2,3, 
for  each  of  the  three  parameter  directions. 

To  improve  efficiency  the  parametric  knot  removal  routine  implemented  is  constructed 
in  a way  that  lets  it  abort  the  computation  if  an  approximation  for  any  component  of  the 
parametric  curve  fails  to  lie  within  the  specified  tolerance.  This  fact  suggests  a search 
strategy  where  we  compute  successive  approximations  to  the  initial  spline  by  adding  one 
interior  knot  at  a time,  starting  with  zero  interior  knots,  and  where  each  intermediate 
approximation  is  given  by  the  first  of  these  approximation  processes  to  be  completed. 
Intuitively  we  would  expect  such  a sequential  search  strategy  to  perform  best  for  “large” 
tolerances  and/or  large  problems,  where  it  is  more  to  gain  by  aborting  an  approximation 
process.  In  this  example  we  have  compared  this  search  strategy  with  a strategy  proposed 
in  [6]  using  a binary  search. 

In  all  the  tests  we  have  used  an  initial  trilinear  spline  constructed  by  sampling  the 
function  given  by  f(x,y,z ) = |[sin(27ra;)  + sin(27n/)  + sin(27T2:)]  in  the  points  specified 
by  a uniform  3-dimensional  grid  on  the  domain  fi  = [0,  l]3,  for  four  selected  grid  sizes. 
Each  spline  was  reduced  by  using  both  of  the  search  strategies  mentioned  above,  for 
tolerances  varying  from  e = 0.001  to  e = 0.01.  Both  of  the  search  strategies  produced 
approximately  the  same  end  grid  size  in  each  test. 

In  Figure  1 the  CPU-time  of  the  two  search  strategies  is  plotted  against  the  tolerance 
for  the  selected  grid  sizes.  We  observe  that  the  reductions  utilizing  a binary  search 
perform  best  on  small  problems,  while  the  sequential  search  strategy  turn  out  to  be 
superior  for  large  problems. 
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(c)  Problem  size  2503 


Fig.  1.  A comparison  of  two  different  search  strategies. 

Example  5.2  In  this  example  we  compare  the  two  different  knot  removal  methods 
presented  in  this  paper.  Here  we  have  used  an  initial  trilinear  spline  constructed  by 
sampling  a function  given  by  f(x,  y,  z)  = cs,n'27rT  in  the  points  specified  by  a uniform 
3-dimensional  grid  on  the  domain  0 = [0,  l]3,  for  varying  grid  sizes.  Each  spline  was 
reduced  by  both  the  method  based  on  the  symmetric  approach  and  the  method  treating 
one  parameter  direction  at  a time. 


The  results  are  presented  in  Table  1.  We  see  that  in  our  implementation  the  method 
using  the  symmetric  approach  is  by  far  the  slowest  method.  However,  at  least  for  the 
type  of  function  considered  in  this  example  the  method  based  on  the  symmetric  approach 
will  give  a much  better  reduction  than  the  other. 
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Knot  Removal  for  Trilinear  Splines,  Tolerance  e = 0.005 

Start 

grid 

Parametric,  binary  search 

Symmetric,  binary  search 

CPU 

End  grid 

Error 

CPU 

End  grid 

Error 

1003 

16.53 

72  x 65  x 65 

4.93800  ■ 10“3 

63.23 

54  x 53  x 53 

4.92080  • 10-3 

1503 

56.44 

81  x 71  x 71 

4.80243  • 10~3 

122.2 

51  X 49  x 49 

4.77236  • 10“3 

2003 

99.48 

68  x 66  x 66 

4.91142-  10“3 

300.9 

54  x 50  x 51 

4.98275  ■ 10-3 

2503 

165.3 

74  x 62  x 62 

4.74970  • 10“3 

584.8 

61  x 56  x 56 

4.85916  ■ 10~3 

300s 

256.8 

72  x 62  x 62 

4.85316  - 10~3 

1094 

60  x 54  x 53 

4.81551  • 10“3 

350s 

391.4 

75  x 65  x 63 

4.77028  ■ 10“3 

1312 

54  x 50  x 50 

4.92422  • 10-3 

4003 

494.6 

71  x 59  x 63 

4.79631  - 10-3 

1865 

54  x 50  x 50 

4.81064  ■ 10~3 

Tab.  1 Knot  removal  for  the  trilinear  splines  of  Example  2. 
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